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ABSTRACT 

We investigate how accretion episodes onto massive black holes power quasars and 
active galactic nuclei while they accumulate mass into the holes. We implement an 
analytic approach to compute both the trend and the stochastic component to the 
trigger of the accretion events, as provided by structure buildup after the hierarchical 
paradigm. We base on host galaxy evolution proceeding from the protogalactic era at 
redshifts z 2.5 dominated by major merging events in high density regions, to the 
subsequent era marked by galaxy-galaxy interactions in newly forming groups. These 
dynamical events perturb the gravitational equilibrium of the gas reservoir in the hosts, 
and trigger recurrent accretion episodes first in the Eddington-limited regime, later in 
a supply-limited mode controlled by energy feedback from the very source emission 
onto the surrounding gas. Depletion of the latter by these events (adding to quiescent 
star formation) concurs with the slowing down of the clustering to cause a fast drop 
of the activity in dense regions. Meanwhile, in the "field" later and rarer events are 
triggered by interactions of still gas-rich galaxies, and eventually by captures of dwarf 
satellite galaxies; these are also included in our analytic model. Thus we compute the 
quasar and AGN luminosity functions; we find these to brighten and rise from z ~ 6 
to z ~ 2.5, and then toward z ~ to dim and fall somewhat, in detailed agreement 
with the observations. From the same accretion history we predict that for z < 2.5 the 
mass distribution of the holes progressively rises and shifts rightwards; we compare our 
results with the local data. We also find that downward of z ~ 2.5 the Eddington ratios 
related to emitting, most massive holes drift below unity on average, with a widening 
scatter; meanwhile, some smaller holes flare up closer to the Eddington limit. We 
conclude that the accretion history, however rich, is dominated by dwindling events 
triggered by interactions under the control by feedback; these establish a link between 
the declining but widely scattered distribution of the Eddington ratios and the tight, 
closely stable upper section of the M — a correlation. Two clearcut predictions arise 
from our interaction picture: the BH activity will appear to follow an anti-hierachical 
trend; the QSO-AGN population is expected to be bimodal, and related to the bimodal 
galaxy population. 

Key words: black hole physics - galaxies: active - galaxies: interactions - galaxies: 
nuclei- quasars: general . 



1 INTRODUCTION 

Accreting black holes (BHs) with masses M ~ 10 6 - 10 9 M e 
are widely held (Rees 1984) to energize Active Galactic Nu- 
clei (AGNs) and quasars (QSOs) and produce their huge 
bolometric outputs that approach L ~ 10 48 erg s" 1 (e.g., 
Groote, Heber & Jordan 1989, Hagen et al. 1992, Maraschi 
& Tavecchio 2003). Not only strong gravity but also large 
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mass inflows M = L / 77c 2 up to some 10 2 Moyr -1 must be 
involved in these objects, even when the overall efficiency 
77 for converting gravitational into radiative energy is up to 
10%. 

By the same token the BH masses 

M ( f ) = -K I dt'L(t') (1) 
VC J 

are expected to keep the archives of the luminous history 
of these sources down to the cosmic time t; barring major 
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silent accretion, the above equation relates the cumulative 
variable M to the quantity L which signals current source 
activity. 

The integral nature of Eq. (1) allows three main activ- 
ity patterns, as originally discussed by Cavaliere & Padovani 
(1988), see also Kembhavi & Narlikar (1999). For one, the 
QSO population could be comprised of a limited number 
of sources continuously active over several Gyrs, that accu- 
mulate large masses M(t) in excess of 10 10 Mq while emit- 
ting lower and lower outputs L(t), to agree with the fall 
observed in the bright QSO population for redshifts z < 2.5 
(see Osmer 2004 for a review). But then the ratios L/M 
ought to drop sharply with the cosmic time t; strong t- 
dependence would also occur in the relation between the 
increasing BH masses and the more persisting properties of 
their host galaxies, e.g., the velocity dispersion a. However, 
this specific pattern is ruled out by the observations of super- 
massive BHs widespread in many local and currently inac- 
tive galaxies, with an upper mass envelope at M ,$ 5 10 9 Mq 
(see Tremaine et al. 2002); a similar bound is also found 
for QSOs shining at higher z (see McLure & Dunlop 2003, 
Vestergaard 2004). 

At the other extreme, Eq. (1) holds as well during a sin- 
gle accretion event of a mass fj, over the time scale r ~ 10" 1 
Gyr to yield a flash of bolometric luminosity L ~ t]c 2 /j,/t. 
Here each BH's mass is accreted in one go over a constant 
r, so M — (i applies and the ratios L/M ~ r\c 2 / r ought 
to be closely constant. Here the QSO fall would be made 
up by many sources active only once and briefly, conceiv- 
ably self-limited by the radiation pressure at the Eddington 
value L — Le = Mc 2 /tE, with M(t) exponentiating on the 
Salpeter timescale rjtE — 5 10 -2 Gyr. Over the cosmic time 
t new BHs ought to form continuously, but with progres- 
sively lower masses for z < 2.5 so as to track the observed 
population decline to lower luminosities while retaining con- 
stant Eddington ratios Xe = L/Le — 1. For this to occur, 
the trend in the hierarchical formation of structures toward 
ever-growing masses would have to be reversed in a closely 
tuned way for the active BHs (see discussion by Merloni 
2004). Then the tight M — a c correlation observed between 
the BH masses M and the bulge velocity dispersions a c (first 
pointed out by Ferrarese & Merritt 2000, Gebhardt et al. 
2000) would evolve markedly, being progressively extended 
toward smaller M for decreasing z, still with a narrow scat- 
ter. 

Finally, the luminous history may involve in many 
bright galaxies a number of recurrent accretion events with 
decreasing outputs; then Eq. (1) yields 



and the BH mass has the form M — ^ fc /i^ with the few first 
contributions dominating. Now the individual BH's masses 
will grow moderately with t, in accord with the intrinsic 
hierarchical trend but in accord also with the decreasing 
average luminosities. Such limited supplies will cause a mild 
drift of the average Eddington ratios below the value Xe — 1, 
but only weak evolution in the M — a c relation. Here we plan 
to discuss the latter two activity patterns, motivated by a 
number of observations. 

Shields et al. (2003) observe the M—a c correlation to be 
largely in place at epochs z — 2.5 at least in its upper range, 



the masses being afterward modified by factors of a few at 
most. The single-flash pattern of activity with Xe — const 
hardly accords with these findings, see discussion by Corbett 
et al. (2003); it also disagrees with the considerable scatter 
shown by ratios A_e. 

As to the latter, Eddington luminosities are welcome 
at the highest observed redshifts z > 6, when the age of 
Universe was definitely shorter than 1 Gyr and yet powerful 
QSOs with L ~ 10 47 erg s" 1 already appear in the SDSS 
data (Fan et al. 2003). These sources require minimal BH 
masses M ~ 1O 9 M0 if they shine at Xe — 1, a condition 
that also fosters fast growth from much smaller seeds. In 
fact, for z — 2.1 McLure & Dunlop (2003) find sources with 
Eddington ratios often attaining values A_b — 1 (see also 
Corbett et al. 2003), but undergoing a slow average decrease 
toward z ~ 0.2. In addition, Vestergaard (2004) finds A_b to 
show large, partly intrinsic scatter at all redshifts z ^ 4, 
with a declining upper envelope for the massive objects as 
z decreases. These observations indicate that the radiation 
pressure limit to accretion is often attained at high z, but 
also that it does not always constitute the tightest constraint 
at lower z. We propose that for z^2a supply limit may be 
more effective for many BHs. 

The counterpart of Eq. (1) extended to the entire un- 
obscured population (or populations) is provided by the re- 
lation (cf. Soltan 1982) 

J dMN(M, t)M = J dt' J dL N(L, t') L . (3) 

Here N(L, i) is the bolometric luminosity function (LF) of 
the quasars and AGNs, and N(M, t) is the mass distribution 
(MD) of the related BHs. 

The development of the former is known to be sharp, 
coherent and non-monotonic; the QSO luminosity density 
peaks at the epoch corresponding to z ~ 2.5 — 3, and rises 
and falls on both sides on timescales of a few Gyrs (Schmidt 
1989, Osmer 2004). In detail, the optical LF "evolves" 
sharply from the Local Universe, with the number of lu- 
minous QSOs rising out to z ~ 2.5 (see Grazian et al. 2000; 
Croom et al. 2004). Past this redshift, "negative evolution" 
occurs with the luminous QSOs dwindling on average out to 
z ~ 5 and beyond (see Kennefick, Djorgowski & De Carvalho 
1995; Schmidt, Schneider & Gunn 1995; Fan et al. 2001). A 
similar trend is also shown by the LF observed in the X-ray 
band out to z ~ 3 (see Miyaji, Hasinger & Schmidt 2000; 
La Franca et al. 2001; Cristiani et al. 2004a). 

Less known is the behaviour of the MD. The local distri- 
bution has been estimated from the observations of a num- 
ber of relic supermassive BHs (see Yu & Tremaine 2002; 
Marconi et al 2004; Shankar et al. 2004). These authors 
find the integrated mass density pth = J dM N(M, t ) M 
to take on local values in the range 4 — 6 1O 5 M0 /Mpc 3 , con- 
sistent after Eq. (3) with the optical emissions from QSOs 
summed to the X-ray emissions from AGNs, with a con- 
siderable contribution from obscured sources. In the range 
2 < 2.5 where the quasar population declines, these AGNs 
still contribute appreciably to the BH masses; eventually, the 
dominant contribution comes from the AGNs pinpointed in 
X-rays (Hasinger 2003, Fabian 2004, Cristiani et al. 2004a). 

We stress that Eqs. (3) and (1) with their integral na- 
ture require neither the variables L and M to be linked 
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by a constant value of Xe, nor the functions N(M,t) and 
N(L, t) to be linked by a relation of the type N(L, i) dL oc 
8 N(M,t) dM in terms of the emission duty cycle 8 ~ r/t. 
In fact, both assumptions, however appealingly simple, run 
at some point against a piece of the available data, espe- 
cially on our side of the QSO peak. In this paper we shall 
relax these assumptions, and will widen our scope to de- 
rive the actual relationship between L and M and between 
N(L,t) and N(M,t); in §2 we give our guidelines. In §3 we 
describe how supermassive BHs form and shine at z > 2.5 
while their protogalactic hosts are built up. In §4 we com- 
pute how formed BHs are re-fueled during the subsequent 
era z < 2.5 by interactions of their hosts with companion 
galaxies in dense regions. In §5 we develop our formalism 
to evolve the BH mass function throughout the two above 
eras. In §6 we derive the QSO luminosity function evolv- 
ing under the drive of the above BH accretion history. In 
§7 we consider the late fueling events due to interactions in 
low density regions and to captures by the host of satellite 
galaxies. In §8 we sum up and discuss our findings. 

Our framework will be provided by the "Concordance 
Cosmology" with H ~ 70 km/s Mpc, fl m ~ 0.27, Qa ~ 
0.73, fib ~ 0.044 (see Bennett et al. 2003); this im- 
plies the relation we use between t and z, namely, t = 
11 sinh' 1 [1.6 (1 + z)" 3/2 ] Gyr. 



2 OUR GUIDELINES 

Many authors have attempted to model the connected evo- 
lution of the BH MD and of the shining QSO LF, particu- 
larly: Cavaliere & Vittorini (2000 and 2002, CV00 and CV02 
hereafter), Kauffmann & Haehnelt (2000), Cattaneo (2001); 
more recently yet, Volonteri, Haardt & Madau (2003), Menci 
et al. (2003), Wyithe & Loeb (2003), Steed & Weinberg 

(2003) , Marconi et al. (2004), Granato et al. (2004), Merloni 

(2004) . 

In the present paper our guideline will be that such a 
connection has to take place through the Eqs. (1) and (3) 
intrinsic to any accretion process. These processes are far 
from trivial, since to power sources as luminous as the bright 
QSOs they have to involve large gaseous masses and high ac- 
cretion rates. For outputs exceeding L ~ 10 47 Lq erg s _1 , 
masses /i ~ 1O 9 M0 or more must be funnelled toward the 
central BH in times r ~ 10~ Gyr or less. In addition, these 
accretion episodes must be coordinated into the highly co- 
herent, non-monotonic evolution the QSOs exhibit on scales 
of a few Gyrs over a span of about 13 Gyr. 

Large gravitational perturbations involving a consider- 
able volume of the host galaxies, such as merging of proto- 
galaxies or galaxy-galaxy interactions, are required in order 
to distort the galactic potential on kpc scales, and remove 
enough angular momentum from the galactic gas so as to 
start it on an inward course toward the central BH. 

But major mergers of the dark matter (DM) host haloes 
as the only accretion triggers hardly can account for the 
observed strong evolution of QSOs for z S; 2.5. In fact, as 
discussed in detail by Menci et al. (2003), the widely agreed 
paradigm of hierarchical growth of cosmic structures implies 
the merging rate at galactic scales to decrease slowly by a 
factor of some 10 _1 from the QSO peak to the local Uni- 
verse. This by itself would lead to overpredict the bright 



QSOs surviving at low z, at variance with the dramatic fall 
observed in the QSO population which entails factors of or- 
der 10~ 2 in the number of bright sources. Even when ac- 
count is taken of the galactic gas exhaustion due both to 
BH accretion itself and to thhe ongoing star formation (see 
Kauffmann & Haehnelt 2000, Cattaneo 2001), the computed 
population of QSOs with blue magnitudes Mb ^ —24 would 
still outnumber that observed for z ^ 1. Moreover, the MD 
accordingly computed differs considerably at z ~ from 
what is inferred on the basis of the Al — a observations and 
the distribution of the velocity dispersions (see Wyithe & 
Loeb 2003). 

We instead relate BH growth with QSO and AGN lumi- 
nosities from considering the rich picture that arises natu- 
rally after the hierarchical formation of cold DM structures; 
this envisages the initially minute density perturbations to 
grow, collapse and virialize under the drive of the gravita- 
tional instability, and in closer detail to merge with similar 
clumps into larger structures (see Peebles 1993). We follow 
this standard course, and consider how the various dynam- 
ical events it predicts affect the dynamical life of the host 
galaxies and trigger different modes of accretion onto the 
central BH. 

These events include the early major mergers that as- 
semble the massive protogalaxies, but also include later and 
milder tidal interactions of the formed hosts with compan- 
ion galaxies in dense environments like groups, or even in 
lower density environments such as the Large Scale Struc- 
tures (LSS); eventually, the hosts cannibalize their retinues 
of satellite galaxies. The interactions in dense environments 
(groups in particular) , while still able to trigger large inflows 
of the galactic gas toward the nucleus, outnumber the bound 
mergers for z < 2.5, consume more gas and decay faster, so 
speeding up the QSO evolution (see Menci et al. 2003). 

Later on for z < 1 "field" processes such as interactions 
of gas-rich hosts in LSS and satellite cannibalism together 
are left as the dominant, if often meager fueling mode, so as 
to be phenomenologically perceived as a later AGN popula- 
tion as opposed to that related to dense environments. 

We let all these dynamical events, with their overall 
trend and their stochastic component, to form or rekindle 
the BHs as they may; we just record the outcomes from the 
integral relationship Eq. (1) of M with L. As discussed in 
detail in §4, our key parameter here will be the gas fraction 
/ destabilized by an interaction, with values around 10% 
based on transfer of angular momentum from the orbital 
motion of the partner to the host gas, and confirmed by 
numerical simulations. 

We also compute N(M, t) and N(L, t) separately but 
consistent with the integral relation Eq. (3), and compare 
these distributions with the observed ones. Such observable 
will be related to the statistics of the interactions based upon 
the mass distributions of galaxies in groups, for which we 
adopt the standard Press & Schechter (1974) mass function. 

We stress recent, direct evidences of activity connected 
with clearly interacting galaxies given by Rifatto et al. 
(2001), Komossa et al. (2003), Ballo et al. (2004), and 
Guainazzi et al. (2005), who report AGNs hosted in both 
galaxies of the interacting systems ESO 202-G23, NGC 6240, 
Arp 299, and ESO509-IG066, respectively. These findings 
complement the extensive, long known body of evidence in- 
dicating that some 30% QSO and strong AGN hosts have 
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close companions or show signs of ongoing interactions, see 
the many single observations referred to in CVOO, and the 
statistics by Bahcall et al. (1997) and by Kauffmann et al. 
(2003). The evidence is even more significant in view of the 
different times conceivably taken by inner fueling and by 
outer disturbances, see the discussions by Beckman (2001) 
and by Tadhunter et al. (2005) , with the references therein. 

One feature specific to the present work is our structural 
inclusion of the effects from the source output itself onto 
the BH; we investigate how such an energy feedback from 
QSO emissions onto the surrounding gas can regulate the 
amount actually accreted and stored into the BHs. In fact, 
accretion is expected to be impaired or utterly halted when 
the gas binding energy in the host is balanced by the energy 
deposited by the source into the surrounding gas during a 
dynamical time, see Silk & Rees (1998). 

Whenever the accretion is so feedback- constrained, the 
balance condition clearly translates into a steep M — a re- 
lation between the mass accreted and the depth of the host 
potential well. In closer detail, we will compute how the 
feedback affects the evolution of the M — a relation, the Ed- 
dington ratios, and the evolving LF and MD. As discussed 
in detail in §3, the key parameter here will be the feedback 
efficiency (f> to deposit source energy into the host gas; we 
will use values <j> ~ v/2c ~ a few % following from momen- 
tum conservation from radiation to gas, and consistent with 
independent lines of data. 



3 FORMING BHS IN PROTOGALACTIC 
HALOES FOR Z > 2.5 

Taking up from CV02, it is convenient to divide the growth 
of the BHs into two main regimes that overlap around z ~ 
2.5, with later additions. 

At early epochs z > 2.5, supermassive BHs grow mainly 
during the major merging events that in dense environ- 
ments build up massive protogalactic haloes of masses Mh ~ 



Mi, 
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10 Mq. By these events large amounts of gas are 



destabilized and funnelled towards the galactic centre to 
be eventually accreted onto the BH. Meanwhile, the same 
events also replenish the host structures with fresh gas sup- 
plies, and so sustain the amounts m of galactic gas at nearly 
cosmic levels m/Mh — fif,/fi m — 0.15. As a result, the ac- 
cretion is often self-limited at nearly Eddington rates with 
L oc M, while protogalactic haloes and BHs grow together 
but not necessarily in a proportional fashion. Thus during 
this era it is appropriate to assume that BHs form with 
masses M = Mi n directly related to the mass Mh of the 
galactic hosting halo. 

We may express Mh in terms of the DM velocity dis- 
persion a = A [GMh I R) 1 ^ 2 , noting that for the standard 
isothermal sphere A ~ 0.7 applies, while for the DM profiles 
by Navarro, Frenk & White (1997) a similar value A ~ 0.6 
holds at the virial radius R. To actually relate Mi n to a, we 
will consider two different models. 

• The unconstrained accretion model (UA) focuses on 
BH coalescence directly following the merging of their host 
haloes; this yields the simple proportionality Mi n ~ 10~ 4 Mh 
(see Haehnelt & Rees 1993, Volonteri, Haardt & Madau 
2003). When expressed in terms of a in units <t„ = 200 km 
s _1 this yields 



31O 7 A/ (^)' 



(4) 



The first relation Min oc a obtains at a fixed virialization 
epoch, considering that a oc M\ ' 3 p 1//6 holds in terms of the 
DM density p in the haloes. The steeper course on the r.h.s. 
obtains from considering (Haehnelt & Kauffmann 2000) that 
the 2-dependence may be approximately rephrased in terms 
of an additional a dependence (plus a wide residual scat- 
ter), because at high z the standard hierarchical scaling 
o~ oc p( n - 1 )/ 6 ( n+3 > ~ p~ 1/2 applies, having used the index 
n ~ —2 for the power spectrum of initial density pertur- 
bations on galactic scales. When translated in terms of of 
the stellar velocity dispersions a c oca 12 " 11 in the galactic 
bulges as indicated by Ferrarese (2002), by Baes et al. (2003) 
and by Pizzella et al. (2004) the resulting average relation 
turns out to be too flat compared with the slope observed 
by Ferrarese & Merritt (2000), Gebhardt et al. (2000), while 
the scatter is still too large. Note that a similar argument ex- 
tended to z < 1, when the Concordance Cosmology modifies 
the scaling laws into a oc p( 1 + n )/ 6 ( n + 3 ) ( would yield instead 
too steep a relation Mi n oc a\. 

• These drawbacks leads us to consider the alternative 
model based on feedback- constrained accretion (FCA) that 
includes the energy balance suggested by Silk & Rees (1998) 
and worked out by Haehnelt, Natarajan & Rees (1998); 
analogous results come from more detailed work, e.g., King 
(2003), Granato et al. (2004), and Lapi, Cavaliere & Menci 
(2005). The balance condition reads 



1 Le td 



GM h 
2R 



(•») 



This applies since gas unbinding and outflow occurs and 
accretion is halted when the fractional quasar output <j> de- 
posited within a dynamical time td ~ R/o~ — 10 8 yr into 
the current gas mass m(t) in the host exceeds the binding 
energy of the gas. We begin with assuming effective val- 
ues 4> — 10 -2 based for radio-quiet QSOs upon momentum 
conservation between radiation and gas that yields (in the 
absence of cooling) values up to v/2c ~ a few %; for radio- 
loud QSOs the kinetic energy in the jets affords higher effi- 
ciencies, but the statistics of such sources is down to 10%, 
so conserves the weighted value. The latter is independently 
confirmed on considering (Lapi, Cavaliere & Menci 2005) its 
effects on the density of, and the X-ray emission from hot 
gas in groups of galaxies surrounding the quasar hosts. We 
also assume cj> to be independent of the host mass, a point 
to be discussed in §8. 

Recall now that the hosts are resupplied with fresh gas 
under halo merging; so as long as the overall gas mass retains 
nearly cosmic values m ~ Mfc,fij,/fi m oc a 3 , the balance 
condition straightforwardly yields 



M lr , 



fib 



fi m 2GA 



if^ 31 o'* e (0 



(6) 



In terms of o c , this is close to the observed slope (the pref- 
actor at z ~ after additional mass has been accreted is 
discussed below and given in Eq. 30). Eq. (6) may be recast 
into the form M m = 6 10~ 6 (1 + z) 5/2 (M h /10 13 M Q ) 2/3 M h , 
to make clear that feedback constrains smaller haloes to 
form or grow BHs with a lower efficiency Mi n /Mh oc 
(M h /10 13 M Q ) 2 / 3 . 
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In both models, during this era the QSO outputs 
- sustained at Eddington levels - are distributed af- 
ter a LF directly related to the halo MD by the re- 
lation N(L,t)dL = Td+N(M h ,t)dM h ; this is in terms 
of the hierarchical halo formation rate df N(Mh,t) oc 
[M c {t)/M h ] {n+3)/3 N(M h ,t)/t, which decreases toward 
masses smaller than the average mass M c (t) virializing at 
t (see Cavaliere, Colafrancesco & Scaramella 1991), so that 
the r.h.s. only approximately reads Nh(Mh,t) dMh r/t. By 
the same token, the hole MD N(Mi n ,t) is also linked as 
given by 



N in (M in ,i) dM in = N h (M h , t) dM h = N h (a, t) da 



(7) 



to Nh{Mh,i), or to the equivalent a distribution N n (o~, i). As 
for the former, we use the simple expression first proposed 
by Press & Schechter (1974) updated to the Concordance 
Cosmology. The MD computed from Eq. (7) constitutes at 
2 ~ 2.5 an initial condition to be evolved afterward under 
the drive of the interactions described in §4. 

The LF at early z had been preliminarily computed and 
discussed by CV02. They stressed that for z > 3 the FCA 
model with the ensuing non-linear stretching L oc Mi n oc 
(l + z) 5 / 2 yields LF shapes generally natter than their 
UA counterparts, and more in tune with the then existing 
data at high z and bright L (Fan et al. 2001). At fainter L 
the present model - with the smaller prefactor correspond- 
ing to Eq. (6), consistent with the additional accretion we 
envisage at later z - yields an even lower LF, in tune with 
the recent observations by Cristiani et al. (2004b), see Fig. 
6. If needed, the LF may be fine-tuned on using the Sheth & 
Tormen (1999) rendition of the halo MD instead of the sim- 
ple Press & Schechter (1974) form, and on adjusting the still 
unsettled amplitude of the initial perturbation spectrum. In 
sum, three main intrinsic features concur to limit the num- 
ber of small BHs active at early z: first, the lower efficiency 
in forming smaller haloes given by d+N(M h ,t), yielded by 
the hierarchical formation; second, the lower efficiency of 
these in forming BHs after Mi„/M n oc M^ 3 , due to the 
feedback process; third, the low prefactor that leaves room 
for the later increase of the BH masses computed next. 



4 REFUELING BHS IN INTERACTING 
HOSTS FOR Z < 2.5 

For z < 2.5 major mergers become rarer and rarer at galactic 
scales; fewer new massive BHs are formed at these epochs so 
their number is conserved to a first approximation (improved 
in §7), but their mass can still grow. Now the prevailing 
dynamical events that trigger accretion are best described 
as interactions between developed galaxies, and these occur 
mainly in the small, dense groups that at these epochs be- 
gin to virialize. By the same token, the gas mass m(z) in 
the hosts is consumed with no fresh imports provided by 
mergers. So the accretion becomes supply-limited and can 
be easily SM&-Eddmgton. 

Small groups with mass exceeding 10 13 Mq , radius Rg 
and bright galaxies membership N g ^ 3 provide particularly 
suitable sites for the hosts to interact with their companions 
(CV00); in fact, in early groups the density n g of galaxies 
is high, while their velocity dispersion V is still comparable 
with the galactic a, conditions that favour effective binary 



interactions. These are mainly in the form of fly-by, that is, 
binary encounters with impact parameter b ranging between 
the galactic radius R and the value Re N g , that need 
not lead to bound mergers. The cross section is still close 
the geometrical value E ~ AitR 2 as long as V ^ 2 a applies 
, and the average time between these events is given by 



TV = l/n 9 S V 



(8) 



the local value is r ro ~ 2 Gyr. As groups merge into rich 
clusters, n g decreases strongly following the density p(z) in 
the host haloes, a trend only partially offset by the limited 
increase of V. As a result, the interaction rate r" 1 (t) declines 
with time following pV oc (1 + z) 1 ' 6 in the Concordance 
Cosmology. 

An interaction of the host with a group companion of 
mass M' will perturb the galactic gravitational potential, 
and destabilize a fraction / of the cold gas mass m in the 
host from its equilibrium at r ~ kpc from the centre. The 
amount fm funneled to the galaxy centre ends up in part 
into circumnuclear starbursts, and in a smaller part trickles 
down to the accretion disk ending up onto the central BH. 
When the main constraint governing the gas equilibrium is 
provided by the angular momentum j, the fraction / may 
be computed as in CV00 to read 



AG 



M' 
bVa' 



(9) 



This ranges from some fmin — 5 10~ 2 (<t»/<t) to fmax — 1/2; 
the latter constitutes the expected maximal gas fraction 
driven into the central 10 2 pc, as confirmed by aimed numer- 
ical simulations of galaxy interactions (see Mihos 1999). Cor- 
responding to larger galaxies being more resistant to gravi- 
tational distortion, Eq. (9) shows / to scale as <r -1 . 

A fraction around ~ 1/10 of the inflowing mass reaches 
the BH rather than ending into circumnuclear starbursts, as 
indicated by the statistics of the energy sources that heat 
up the dust in bright IR galaxies (see Franceschini, Braito 
& Fadda 2003). So p = fm / 10 is the mass made available 
for actual accretion, while the rest ends up into stars or is 
dispersed. The process of fueling takes times of order td, the 
host dynamical time, and spans a few Salpeter times. 

Between the limits fmin , fmax the probability density 
for / due to the distribution of the orbital parameters pri- 
marily reflects the distribution of the masses M' of the in- 
teraction partners; this is because the encounter velocity V 
and the impact parameter b vary in narrow ranges, and are 
actually correlated in a galaxy group. So / oc M' closely 
applies; since M' is distributed after p(M') oc M'~ 3 with 
s slightly under 2, following the Press & Schechter (1974) 
distribution in its power-law section, the result is close to 
p(f) oc f~ 2 . Th result reads 

p(f\o) = /"°*/"»" /-^ for fmm <; f ^ fmax . (10 ) 

Jmax Jmin 

This is used to compute the average value (/} ~ 15%; since 
/ oc fi applies at given m, the result closely reads p{n) oc 

Recurrent interactions will iteratively exhaust the initial 
gas mass mm in the host; after q interactions the residual 
mass reads 
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m in - f k ) 



(11) 



where each fk is extracted from the probability distribution 
Eq. (10). 

We will often use the simple estimate for the average 
value (m q ) ~ m in (1 — (/))*, based on equal average deple- 
tion factors (/}; this is accurate to within 0.05 nii„ up to 
q — 7 steps as we shall find below to occur on average from 
z ~ 2.5 to 0.2. We may also write for the average depletion 
rate the equivalent differential equation 



rh/m ~ - </)/r r 



(12) 



In fact, its discretized solution starting from the initial con- 
dition rrii n reads just as (m q ) above when computed at the 
step q given by the integer part of 



dt/Tr{t)\ , 



(13) 



starting from tin that corresponds to z = 2.5; the r.m.s. 
deviation of the stochastic variable q is about y/q to a suffi- 
cient approximation. To the depletion given by Eqs. (11) or 
(12) we systematically add the gas consumption by ongoing 
star formation, that we compute following Guiderdoni et al. 
(1998). 

As for the mass actually accreted, in this regime of 
supply-limited accretion we again consider two possibilities. 

• The unconstrained accretion model (UA), where the 
accretion is not affected by dynamical feedback, but still are 
subject to the radiation pressure limit. The mass accreted 
in each interaction is p = fm(z) / 10, distributed as given 
by the simple counterpart of Eq. (10), namely 

p{n\cr) = g M~ 2 f or fmin m/10 < /4 < fmaxTn/10, (14) 

where g = 10 _1 m f max fmin /(,/ 'max — fmin) provides the 
normalization. 

• The feedback- constrained accretion model (FCA), in 
which a tighter constraint to pi is set by the analogous of Eq. 
(5). This is evaluated on considering that now the accretion 
may be sub-Eddington, so Le is to be replaced by r\ fic 2 / r, 
to read n/m ~ (cr/c) 2 r / 2 A 2 td V 4>\ so the constraint reads 
now fj, ^ fxi = m (a/c) 2 (rj (f)) -1 , considering that r ~ td 
and A 2 ~ 0.5 apply. 

The mass actually accreted is given by the minimum 
between the amount /j, = /m/10 made available by the 
interaction, and the constraint set by the feedback; that is 
to say, 



H — min [m //10 , fj,i] 



(15) 



After q interactions m is rescaled down iteratively as said, 
to yield the estimate pn = Mi n (l — {f)) q r/ntE- Note that 
the total probability of the value /i; is contributed by all in- 
teractions leading to accretion events that, if unconstrained, 
would exceed this value; this leads to piling up of accretion 
episodes at the upper bound /i;. To account for this, we write 
the counterpart of Eq. (14) in the form 

p(fJ.\cr) = 5(/i— /i;) for Wfit^fminm 

= g[/i~ 2 + g S(/i — in)] otherwise (16) 
where g' = {fmax m— 10^;) / (f ma x m fit). 

In all cases the luminosity L oc /i attained in any one 



accretion event no longer is in a fixed relation to the current 
BH mass M. This is because the accreted mass fj, depends 
now on stochastic orbital parameters as given in Eq. (9), 
and is distributed according to Eq. (14) or Eq. (16) in the 
UA or the FCA model, respectively. 

We end this Section by giving a simple estimate of the 
final mass of a BH after a number q of interactions; when 
the feedback constraint given by Eq. (15) is effective so as 
to join smoothly at z = 2.5 with the mass Mi„ similarly 
constrained by Eq. (6), a simple upper bound is given by 



M in + M in £(!-</»* 



= M in [l-(l-</»* +1 ]/</> . (17) 
In fact, this is close to the actual value in small and inter- 
mediate haloes when the feedback constrains fi to be close 
to m, see Eq. (15). 

With (/} = 15% we find At ^ 6M in , a bound actu- 
ally approached when the interaction number q grows large. 
More realistically, as q is related to t (or z) by Eq. (13) on 
average, the host undergoes q ~ 7 interactions from z = 2.5 
to 0.2, of which only the initial 4 or 5 are effective on aver- 
age; correspondingly, the masses grow by a factor up to 4. 
On the other hand, the mass remains unchanged for the BHs 
that were never re-activated after z = 2.5; so the growth of 
M/Min spans the range from 1 to 4 at the outmost. This 
may be rephrased in terms of an overall scatter bounded by 
a factor 2 from the average, that is, log M is bounded by 0.3 
dex. 



5 EVOLVING THE BHS FOR Z < 2.5 

On long time scales t > r r the development of the MD 
of the BHs may be viewed at as a stochastic process that 
increases M, at given a and given halo distribution Nh{cr)\ 
the corresponding distribution is denoted by N(M, o~,t). 
This is ruled by the equation 

d t N(M,a,t) = -— N(M,a,t) + 

Tr 

+ — dfip(fi\a)N(M - n,a,t) . (18) 
r r J 

proposed by CV02, and used also by Yu & Tremaine (2002), 

Hosokawa (2004), Menou & Haiman (2004). The evolution- 
ary rate dtN is contributed by two terms. The first describes 
the BHs which interact and thereby increase their initial 
mass M, so depleting the number N(M, a) dM in the mass 
range (M — M + dM). The second describes the number 
of BHs which start from a lower mass M — pi and accrete a 
gas amount /i, with probability p(fi, o) given by Eq. (14) or 
(16) for the UA or r the FCA model, respectively. Here r r is 
the average time between two subsequent interactions of a 
galaxy, discussed in §4; moreover, a w 0.3 is is the host frac- 
tion in dense environment, corresponding to the 30% bright 
galaxies residing in groups with membership > 3 (Ramella 
et al. 1999). In the above equation the number of BHs is con- 
served, while they are re-distributed toward larger masses; 
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to a next approximation number conservation may be re- 
laxed on adding to Eq. (18) the appropriate source term as 
discussed in §7. 

To capture the evolutionary trends given by Eq. (18) it 
is convenient to consider at first small accretion events with 
/x/M << 1, and Taylor expand to second order. So, we end 
up with the approximate equation 

d t N{M,a,t) ~ _^W_ d M N(M,cr,t) 

+ d 2 M N(M,a,t). (19) 

This is similar to a Fokker-Planck equation, actually one 

based on the probability distribution p(fi j a). The coefficient 
of the first order derivative C(Al,t) = (/i) / r r represents 
the average upward drift of the mass under accretion, while 
the coefficient of the second derivative D(M, t) = {/j, 2 } / 2r r 
plays the role of a diffusion coefficient; the averages are com- 
puted on using the probability distribution p([i \ a) given in 
Eq. (14) or (16). Note that in the context of QSO evolu- 
tion (Cavaliere et al. 1983, Small & Blandford 1992), Eq. 
(19) constitutes a continuity equation for the MD that con- 
tains also a diffusive term; correspondingly, N{M, t) not only 
drifts toward larger masses, but it is also reshaped reflect- 
ing the scatter in the masses added by the stochastic re- 
activations. 

Here we give the analytic solution of the above equa- 
tion in the simple case when the coefficient D(M,t) is con- 
stant in time and independent of M (see Shankar 2001). 
We solve for the evolution from an initial mass distribu- 
tion aNi„(M,a) of the BHs residing in groups at z ~ 2.5. 
We perform a transformation to Lagrangean coordinates 
M c = M — J C(M,x)dx, where the drift term is absorbed 
into a total derivative; so we end up with a pure diffusive 
equation. This is solved by standard methods to yield 

N{M c ,a,t)= ° / d£Ni n (£,a)e lM M* . (20) 
2v %Dt J 

Finally, we go back to Eulerian coordinates and obtain 
our solution N(M,a,t). We represent in Fig. 3 the local 
mass distribution N(M,t ) = J da[(l - a)N in (M,a) + 
N(M,a,to)], with the first term due to the dormant BHs 
that do not reside in groups; both terms are integrated over 
the variable a. We use the following parameter values: the 
coefficient C(M,t) ~ 10 s (Af/10 9 M ) 2/3 (t/t )~ 2 M /Gyr 
is derived from its definition combined with Eqs. (9) and 
(12), averaged with the use of Eq. (10), and considering 
also that m oc Mh oc er 3 applies; the second coefficient 
D ~ 310 14 M Q /Gyr is derived similarly, but performing a 
final time average. 

The solution of Eq. (19) in the form M N(M,t) is rep- 
resented by the solid, thin line in Fig. 4. It is seen that the 
evolutionary trends described by the approximate solution 
are close to those from the full Eq. (18) where the actual 
mass additions fj./M are finite; the similarity is closer in the 
FCA case where the probability in Eq. (16) is effectively 
confined to a narrow range. 

To solve the full Eq. (18), we use recursive stepwise inte- 
gration, having assumed all host galaxies residing in groups 
(a fraction a of the hosts) to undergo an interaction in the 



time interval t t . The BHs are formed with initial mass Mm 
by the end of the era z 2.5; each of the subsequent interac- 
tions (labeled by the index k) contributes an additional mass 
pu- After q interactions the BH mass is given by the sum 
M — X^fc-o ^ k °f the stochastic amounts /Ltfc, with /io = Mi„. 
The probability density for p^ is given by Eq. (14) for uncon- 
strained accretion; if instead the constraint by the feedback 
is effective, [ik will have the probability density given by Eq. 
(16). 

Thus after q interactions the solution for N(M, a, t) at 
the time t provided by Eq. (13) is given by 

N(M,a,t) = N h {a) [(1 - a)P ln {M \ a) + aP q {M \ a)] . (21) 

Here again the first contribution is due to the dormant BHs; 
P q (M\a), the conditional probability to find a BH mass 
M in a halo with given a, is computed with the recursive 
equation 

P q (M \a) = J dp, p(fj, | a) P q -i(M - ft | a) . (22) 

This starts out with the conditional probability for the initial 
step k = 0, that we express as 

Pin(M | a) = S[M - M in {a)} (23) 

by continuity with the M — a correlation produced at the 
end of the previous era z > 2.5. 

Each sheet in Figs, la and 2a represents the contribu- 
tion, expressed by Eq. (21), to the MD of relic BHs from 
hosts with a given velocity dispersion a. The contribution 
from hosts that never previously interacted is represented in 
the form of spikes peaked around the BH masses Mi n formed 
at z ^ 2.5. The alignment of such spikes shows how these 
masses are related to a after the Eqs. (4) or (6), respec- 
tively. At ~ 2.5 the regime of galaxy interactions in groups 
begins; now the spikes partially drift and spread along the 
mass scale, as the BHs grow by such stochastic accretion 
events. The growth differs in the UA and FCA model. 

In the unconstrained model UA (Fig. la) the BHs grow 
by stochastic amounts fj, distributed with the probability 
p{n \ o) given by Eq. (14). This extends over the mass axis 
with decreasing values. In the feedback- constrained model 
FCA (Fig. 2a) the growth of the BHs is reduced, being 
bounded by a factor 4. This is because the constraint cuts 
off the upper range of the probability distribution p(/x, a) 
tucking it - as it were - at the upper bound of the range of 
(i, as described in detail by Eq. (16). 

In Figs, lb and 2b we project onto the M, a plane 
the probability N(M, a,z = 0)/N h (a) given by Eq. (21); we 
have also reported the data concerning M—a as discussed by 
Tremaine et al. (2002) . We predict the feedback-constrained 
mass accretion to be not so abundant as to materially change 
the early M — a correlation. So in our constrained model 
this is rooted back at z ^ 2.5, with only mild alterations 
occurring afterwards, and these mainly at intermediate and 
small masses; the result is consistent with the observations 
by Shields et al. (2003). 

In Figs. 3 and 4 we show the behaviour of N(M, t) — 
J daN(M,a,z), the MD integrated over o, i.e., over the 
sheets shown in the previous two Figures; note the depen- 
dencies on a of Eq. (21). We illustrate how the MD changes 
from z — 6 (dashed lines) to z = 2.5 (dotted lines), and then 
to z — (thick solid lines) ; in both Figures the shaded region 
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represents the observational evaluations from local data in 
different bands (Yu & Tremaine 2002, Shankar et al. 2004). 
In the UA model shown in Fig. 3 the relic MD at z ~ ex- 
tends toward large masses with a high tail; here the match to 
the data is poor, since the latter indicate a sharper decline. 
Fig. 4 shows how in the FCA model the feedback depletes 
the tail and produces a closer fit to the local data; the accre- 
tion so constrained goes to contributing significantly more 
in the intermediate mass range where it causes the MD to 
swell, again in agreement with the data. We have included 
the contribution of later accretion events occurring in the 
field (to be discussed in §7), to the effect of yet improving 
at small masses the agreement with the observations, as can 
be seen by comparison with the dashed-dotted line which 
does not include this contribution. The approximate diffu- 
sive solution for the MD is represented by the thin solid line 
in Fig. 4. 

The same interactions that produce the MD above when 
time-integrated over several Gyrs, also produce over times 
10" 1 Gyr the LF that we discuss next. 



6 THE LUMINOUS EVOLUTION FOR Z < 2.5 

At high redshifts z > 3 where the luminosities are expected 
to be Eddington limited and feedback constrained as dis- 
cussed in §3, the LF has been computed by CV02; we just 
recall in Fig. 6 their prediction at z ~ 4.5 to show the agree- 
ment with the recent observations by the GOODS survey 
(Cristiani et al. 2004b). 

Here we focus on the the range z < 2.5, an era when 
the supply is limited due to progressive exhaustion of the 
gas reservoirs in the hosts by the many accretion episodes 
that produce L. Two reasons concur to cause here a com- 
plex relation between LF and MD. First, as anticipated in 
§1 an accretion episode of a mass fi over a time r produces 
the luminous output L given by Eq. (2); instead, the cu- 
mulative mass growth is given by M = X^fc^fci involving 
time-integration over the history of the QSO population. 
Second, the mass p, made available for the accretion is ruled 
by the interactions with their stochastic component; so [i 
may be quite smaller than M, and the accretions may be 
very sub-Eddington. In other words, the relation of L tx [i 
with M is no longer tight. 

The actual relation stems from the conditional probabil- 
ity distribution p(p \ a) for accreting \x at given a, expressed 
in Eq. (14) or (16) for the UA or the FCA model, respec- 
tively. This yields also the conditional probability density 
p(L | cr) for the luminosity to attain in a reactivation the 
value L given by Eq. (2), and then fade out. So at any given 
time we have 

p(L|a,t)=p(H<r)^; (24) 

where p tx m{t) effectively depends on t due to Eq. (12). 

Note that in the UA model, the mass M accumulated at 
time t into a BH obeys the relation M(t) = [rrn n — m(£)]/10 
in terms of the residual gas mass m(t) in the host; so, the 
distribution of the latter is directly linked to the MD of BHs. 
This no longer holds in the FCA model, where we use the 
average value for m{i) given by Eq. (12). 

The rate of the reactivations is f3Nf,h/T r , in terms of 



BH number distribution N bh (a,t) = J dM N(M,a,t), and 
of /3 = a/Ng w a/3 that represents the fraction of hosts 
residing in groups (see §5) and interacting over a mean time 
r r ~ a few Gyrs. 

The above components can be brought together to yield 
the LF, upon using the formalism of the continuity equation 
along the L axis, as developed by Cavaliere et al. (1983). 
This takes on the form 

d t N(L, a,t) + L d L N{L, *,t) = f> ^MM^*) , (25) 

Tr(t) 

considering for simplicity light curves equal and monoton- 
ically decreasing on the scale r. The solution is given by 

N(L, a, t) = i jf" dL> p { f) N^WW'WJ) . (26) 

for the numerical computations represented in Figs. 5 and 

6 we use the specific values L — —L/t with r ~ td = 10 _1 
Gyr, while j3 saturates to the value 0.1 soon below z = 2.5. 

Dependencies on a arise in Eq. (26) from the probabil- 
ity distribution p(fi \ a) that enters Eq. (24), and from the 
integrated N\,h((r, t) obtained from Eq. (21). Upon convolv- 
ing over a, we obtain the bolometric LF; this is converted 
to optical luminosity Lb on using the standard bolometric 
correction of 10. The results are plotted in Figs. 5 and 6 
for unconstrained and feedback-constrained accretion - our 
model UA and FCA, respectively - and are discussed below; 
we compare them with the data of Boyle et al. (2000) and 
Grazian et al. (2000). As to z > 3, in Fig. 6 and its caption 
we also recall that in our FCA model the LF has flat shape, 
low normalization, and z behaviour in detailed agreement 
with the observations over a wide range of L at z ~ 4 — 5. 

Points to be noted are as follows. First, the shape of 
the LF may be schematically rendered as a rather flat faint 
section going over to a steeper bright section. The faint sec- 
tion results from the flow along the L axis (2 nd term on 
l.h.s. of Eq. 25) of sources fading from their top luminosi- 
ties (stochastically provided after Eq. 24); the bright section 
results from the shape of p(L) convolved over a. 

Second, the evolution we expect for the LF may be un- 
derstood as a combination of a mild "density evolution" and 
a stronger "luminosity evolution". The former arises as the 
interactions re-activate the BHs in groups on a time scale 
r r that grows moderately with cosmic time; this causes a 
decrease of the amplitude of the LF by a factor around 2 
down to z ~ 0.5. The latter evolution occurs because on the 
same time scale the interactions exhaust the gas content in 
the hosts according to Eq. (12); less remaining gas means 
lower luminosities following L oc m{z). 

Third, the accretion history we envisage produces at 
late z an excess of sources compared to the optically se- 
lected AGNs; this leaves room for sources in X-rays and for 
obscured objects (Fiore et al. 2003, Ueda et al. 2003). A 
larger excess is produced at fainter L by the field processes 
to be discussed in §7. 

Decreasing luminosities and increasing masses produce 
Eddington ratios Xe declining on average. In fact, the lu- 
minosities of the sources when re-activated at z ~ 0.2 are 
lower than at z ~ 2.5 by a factor (1 — {f)) q — 1/3 for q — 7; 
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meanwhile, the masses grow at most by a factor 4 as we 
have seen at the end of §4, so the Eddington ratios decrease 
to about 1/12 toward z ~ 0.2. The actual distribution of 
\e at given a is computed on convolving the probability for 
a luminosity L given by Eq. (24) and that for a mass M 
given by Eq. (22), with the variables combined as to yield 
Xe- This may be represented as 

n 9 (A B ,a) = (27) 
J dM J d f j,P t (M\a)p( f j,\<T)5(X E -n/M), 

considering that for any value of the Eddington ratio fi/M = 
LtE/Mc 2 holds. Fig. 7 illustrates our numerical results with 
feedback-constrained accretion, the FCA model that enters 
through the expressions of p from Eq. (16) and P from Eq. 
(22). 

The figure shows the complex behaviour of A_b that we 
find for increasing number of interactions undergone by the 
hosts. In fact, we have simplified our plots starting them 
from a sharp initial condition at z = 2.5. The overall decline 
is due to the average dimming of all luminosities, more rapid 
in hosts with smaller a; these undergo relatively faster gas 
consumption due to the scaling (/) oc cr _1 discussed in §4. 
The overall range of / is affected by the varying lower end 
fmin oc a -1 and also by the distribution of M'\ the result 
is illustrated by the pairs of lines corresponding to given 
values of a, with two instances named in the caption. The 
outcome is a correlation between L and M which is tighter 
in smaller hosts, more sensitive (in the absence of the lim- 
its from cooling or optical depth as discussed in §8) to the 
feedback constraint. 

In sum, as z decreases from 2.5 the average values of 
(Xe) decline for the truly supermassive BHs in old spheroids 
while the intrinsic range covered by Xe widens, consistent 
with the observations by McLure & Dunlop (2003) and by 
Vestergaard (2004). Next we consider the higher values of 
Xe contributed by interactions involving hosts located in the 
field. 



7 LATER ACCRETION EVENTS IN THE 
FIELD 

To now we have focused on accretion events driven for 
z < 2.5 by interactions of the old host galaxies residing in 
dense environments like the groups. But hierarchical clus- 
tering also envisages additional ways to feed BHs, which for 
2 ;$ 0.5 contribute to accretion onto BHs even when their 
hosts are located in lower density environments (the "field"). 

There some new BHs continue to form by major merg- 
ing events even at z < 0.5, though at a reduced rate; mean- 
while, starving BHs are still refueled by interactions, though 
on longer time scales r/ ~ 10 Gyr. By the same token, these 
later and rarer events involve hosts that are still moderately 
rich of gas, the latter having been depleted only or mainly by 
quiescent star formation; on similar scales, satellite galaxies 
begin their plunge into the hosts, importing some fresh gas. 

So these later, generally smaller but still considerable 
accretion events feed what is perceived as a later population 
of AGNs peaking under z ~ 1. Thereby the LF is enhanced 
especially at faint bolometric L, and is best observed in the 
X-ray band; the integrated contribution to the BH masses 



is appreciable. Since our basic equations Eqs. (18) and (25) 
are linear, in our computations of the observables we have 
summed the outcomes of the following three independent 
processes. 

i) Some major galaxy mergers still occur for z < 2.5, 
if at lower and lower rates. These strong dynamical events 
may form/grow BHs in recently reshuffled, gas-rich hosts at 
low z. Correspondingly, we add a source term to the r.h.s. 
of Eq. (18), in the form P in (M\a) d + N h (a, t) that yields the 
host merging rate in terms of a. This yields some 20% of 
all AGNs at z ~ 0.5; but these "new" BHs will start below 
the overall M — a relation, and shine close to A_b — 1 at 
intermediate Lb- 

ii) Interactions occur also in non-virialized LSS for 
z ;$ 0.5. Here the relative velocities V are up to 300 km s -1 , 
the upper bound to the pairwise ones; the galaxy densities n g 
are lower by about 3 10 -2 , and so are the related encounter 
rates tJ 1 ~ n 9 EV. But the hosts so involved are about 4 
times more numerous than in groups, which include only a 
fraction a ~ 0.3 of the galaxies in the field (Ramella et al. 
1999). The corresponding rate in LSS to be inserted on the 
r.h.s. of Eq. (25) is given by (1- a)N bh (a,t) p(L\ a,t) /r f (t); 
at z ~ 0.5 this contributes about 10% compared with the 
rate in groups, a fraction increasing to 30% toward z ~ at 
Lb ~ 10 45 erg s _1 . On the other hand, in these conditions 
the gas exhaustion rate ra/m ~ (/)/t/ is lower by 3 10 -2 
compared with hosts in groups, see Eq. (12); this leaves in 
the hosts correspondingly more residual gas, thus the QSOs 
so reactivated may burst out at higher Eddington ratios than 
their coeval counterparts in groups. 

iii) The DM merging history also includes many events 
where the hosts end up cannibalizing their satellite galaxies 
(see Menci et al. 2003) together with the associated, scant 
gaseous content; many traces of ongoing such events are be- 
ing unveiled in our Local Group, see Martin et al. (2004); 
Law, Johnston & Majewski (2005). Such episodes cause only 
small accretion and weak, often sub-Eddington AGN emis- 
sions; these are easily drowned into the starlight or obscured 
by dust in the optical band, but are noticeable in X-rays (see 
Di Matteo et al. 2000). On the other hand, the capture rate 
is considerable just under z ~ 1, as the cosmic time ap- 
proaches the time scale of dynamical friction which sets the 
beginning of the satellite plunge into the central galaxy. 

In detail, we consider that the satellites involved have 
masses M s ~ 10 s — 10 10 Mq, and begin their infall into 
the host potential well under dynamical friction in standard 
times t s ~ 3 (M» / 10 10 M )~°' 7 Gyr, see Binney & Tremaine 
(1987). But during the capture the gas in the satellite is 
peeled off by tidal disruption and stripping (Colpi, Mayer 
& Governato 1999); so the gas masses reaching the host 
nuclear region are reduced down to p, s ~ 10 _2 M S . Assuming 
that the initial mass function of the satellites follows the 
low-mass end of the Press & Schechter expression close to 
N S (M S ) oc Mj" 2 , the probability distribution for p, s is given 
by ps(p-s) = 10 B M Q fij 2 in the range 10 6 < fi s < 10 8 M©. 

Thus satellite captures constitute another, late stochas- 
tic process that contributes to accretion. They affect the 
evolution of the MD in a way still described by an equation 
like Eq. (18); but now the average increment p is roughly 
constant in time rather than being proportional to m(t), and 
occurs in many galaxies. Meanwhile, the number of satellites 
is depleted from an initial value around 10, at the rate 
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N(Ms)/N(M s ) = -1/t, , (28) 

so that N(M S ) decreases on scales of several Gyrs. 

We illustrate in Fig. 6 how (ii) and (iii) contribute about 
70% to the faint end of the integral LF at z ~ 0.5, and up 
to 80% at 2 ~ 0. The latter, dominant contribution under- 
goes density evolution on the scales of several Gyrs by the 
decrease of N(M S ); this will be luminosity-dependent owing 
to the basic dependence r s oc M~ ' 7 on the satellite mass. 
Meanwhile, the contribution to the integrated masses from 
(ii) and (iii) is considerable, as shown in Fig. 4. 

Finally, low outputs are conceivably contributed by in- 
dependent processes. In elliptical hosts the cooling of the 
galactic gas can provide sufficient accretion to power low- 
luminosity radio sources, see Best et al. (2005). In spiral 
hosts disc instabilities and bar formation provide enough gas 
inflow to power faint optical AGNs (see Sellwood & Moore 
1999, Combes 2003); it will be interesting to see what such 
processes contribute to the overall LF and MD. 



8 DISCUSSION AND CONCLUSIONS 

We have discussed how the standard hierarchical paradigm 
for structure formation gives rise to a unified if rich picture 
for the accretion history of the supermassive BHs energiz- 
ing the QSO and AGN emissions. The paradigm envisages 
the host galaxies to undergo a sequence of dynamical events 
beginning with early major mergers in high density environ- 
ments, and passing over by z ~ 2.5 to milder interactions 
with companion galaxies in groups; later on, interactions of 
hosts the field join in, and the overall sequence ends with 
captures from the retinues of satellite galaxies. 

All such dynamical events trigger some gas inflow to- 
ward the nucleus, albeit along the sequence the accreted 
masses decrease on average. At the two extremes of the se- 
quence, major mergers and satellite captures both import 
fresh gas into the hosts; while the intermediate interactions 
just tap diminishing gas reservoirs by perturbing the sym- 
metry of the host gravitational potential and inducing non- 
conservation of the integrals that control the gas equilib- 
rium, such as the angular momentum. 

On this basis we have computed several linked observ- 
ables on using a formalism focused on following both the 
trend and the stochastic components to the accretion events. 
We have computed how the statistics of merging events and 
interactions produce in nearly real time (i.e., within lO^ 1 
Gyr) the shape of the QSO luminosity functions N(L). We 
have also computed how their petering out over scales of 
a few Gyrs concurs with the exhaustion of the galactic gas 
reservoirs to produce strong evolution in the LF described by 
our N(L, z) in Fig. 6. Over yet longer scales of several Gyrs 
their cumulative actions add up to grow the BH masses M 
and to change their distribution N(M, z) as shown by Fig. 
4. 

We have seen that in the early era z <; 2.5, when large 
proto-spheroids were in the process of buildup through ma- 
jor merging events, the BH masses closely tracked those of 
the host haloes; correspondingly, the MD of the BHs also 
tracked closely that of the host DM haloes, see Eqs. (7). 
These conditions hold at high z, when not only the proto- 
galactic dynamical scales are close to the Salpeter time, but 
also replenishment of galactic gas is granted by the same 



merging events that trigger the accretion; then the Edding- 
ton limit applies to yield L oc M. In turn, the MD and the LF 
grow together at these high redshifts, to yield approximately 
N(L, z) dL ~ N(M, z) dM r/t, as envisaged by Marconi et 
al. (2004) improving on Small & Blandford (1992). 

Our main focus was the later era z ^ 2.5, when these 
simple relations break down, and L takes on a different 
course from M as does the LF from the MD; basically, this 
is because now a supply limit applies. In fact, the hierar- 
chical paradigm indicates as main triggers the interactions 
of the host galaxies within the small, dense groups that at 
these epochs begin to virialize. These interactions, while no 
longer providing fresh gas to the host, can funnel toward the 
nucleus fractions [i/m of the residual gas m(z) left over by 
previous events, that are considerable yet lower than the BH 
mass already accumulated. So in this era it is convenient to 
represent the luminosities as 

L oc / m(z) , (29) 

to highlight the historic trend embodied in the residual gas 
mass m(z), and the stochastic component / = /i/m trig- 
gered by the last accretion event. At given m the strength 
of such events and the levels of the associated L are set by 
the stochastic distribution p(fi) related to the orbital pa- 
rameters of the interactions. The feedback constraint, when 
effective, cuts off the upper range of /i and tucks in - as it 
were - the distribution p(fi) at its upper end as shown by 
Fig. 2a. 

The result is a basic shape N(L) L oc p(/i) oc L~ 2 mod- 
ified by convolution over a, see Fig. 6. For z < 2.5 the 
z-depending LF embarks on a fast decrease; it is intrinsic 
to our view that the peak of the QSO evolution should be 
found at z ~ 2.5 close to the beginning of the virialization 
era for small and dense groups, the sites most conducive to 
interactions. 

We stress why physical continuity between the two main 
regimes of accretion is bound to arise across z = 2.5, and 
how this is implemented in our analytic calculations, in par- 
ticular as for N(L,z) and N(M,z). The first issue clearly 
goes back to the smooth transition large galaxies — > small 
groups in the hierarchical scenario, that in particular implies 
r r (2.5) ~ td yn (2.5). The second issue is related in the UC 
model to the close equality L oc fii ^ Mj ra , that is ensured 
when the mass actually accreted on the first interaction sat- 
isfies //10 ~ Mi n /m; this is the case with the average values 
/ w 0.1 that we independently evaluate from Eqs. (9) and 
(10). The related semi-analytic computations by Menci et al. 
2003 visualize such a smooth transition. On the other hand, 
in our FCA model the coupling efficiency <f> is clearly contin- 
uous across z = 2.5, and this ensures the condition /ii ^ Mi n 
to hold, as shown in Sect. 4, by Eq. (17) and preceding lines. 
Note that our coupling levels <f) ~ 10~ 2 (independently mo- 
tivated at the end of §2), imply a supermassive BH (active of 
dormant) to inhabit nearly all bright galaxies, with masses 
satisfying the M — a relation. Moreover, Lapi et al. 2005 
show that such values for <j>, when used in the feedback bal- 
ance extended to the hot gas pervading groups, also yields 
the appropriate Lx — T relation for the associated X-ray 
emissions. 

While the first population of QSOs and powerful AGNs 
of high density ancestry is on its decline, at z ~ 0.5 our 
LF shows an excess of faint QSOs over standard Type 1 
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sources; as shown in Fig. 6, the excess relatively increases at 
fainter bolometric luminosities L ^ 10 44 erg s _1 and toward 
lower z, mainly due to the field processes discussed in §7. 
These summed excesses emerge at z ^ 1 as a later population 
of AGNs, consistent with the observed numbers of X-ray 
selected and obscured AGNs, see Hasinger (2003), Fabian 
(2004) 1 

Correspondingly, albeit in a fashion softened by its 
time-integrated nature, N(M, z) changes for z < 2.5 by 
drifting to larger M while undergoing considerable reshap- 
ing. As shown by Fig. 4 the latter is in the form of swelling 
in the intermediate range, related to the stochastic distribu- 
tion p(fi) being tucked in by the feedback limit. The overall 
change we compute ends up in a shape agreeing with the 
local observations; Fig. 4 shows this to be the case through- 
out the observed range, and predicts what will be observed 
at higher z. 

The relic, local mass density in BHs is computed on the 
basis of Eq. (3) from our feedback-limited MD increases, 
and yields a factor close to 3 from z = 2.5 to 0, to attain a 
local value of pbh — 5 lO 5 A/0 Mpc -3 in agreement with the 
estimates by Tremaine et al. et al. (2002) and by Shankar et 
al. (2004) . This includes the excess we find at intermediate L 
and z over the Type 1 AGNs; the value goes up to 6 10 5 Mq 
Mpc -3 on including also the late accretion processes in the 
field. 

Two more imprints are left by this rich history. One is 
to be found in the steep and tight relation of M to the DM 
velocity dispersion a in its upper section. In the FCA model 
we find that Eq. (6) is updated to z ~ by the additional 
mass increase given by Eq. (17), to read 

M = 1.31O 8 M , (30) 

illustrated by Fig. 2b. This translates into M oc a 4 2 - 4 ' 5 in 
terms of the velocity dispersion a c of stars in the bulge (see 
§3), which agrees with the observations (recall that ct, = 200 
km s _1 ). Here the scatter is moderate, in fact, bounded by 
0.3 dex, first because the next recurrent events contribute 
progressively less mass on average; second, because the ac- 
tual accretion is constrained by the feedback. In fact, the 
latter effect is dominant in hosts with a ^$ 300 km s" 1 when 
the quasar output is coupled to the surrounding medium at 
levels <f) ss 10~ 2 . In full, the limit is provided by Eqs. (15) 
to read (er/c) 2 < (5 10~ 2 r\ 4> fmin t d /r); on recalling from §4 
and §3 that fmin = 510 _2 cr*/<7 applies, the condition may 
be recast to read 



a < 300 



Tj 4> fmin td' 


1/3 


r A i 


10- 1 10- 2 0.05 r_ 




_0.6j 



At the other extreme of small haloes we note that a BH 
may end up its trajectory on the M — a plane above the 
values given by Eq. (30) if the coupling is weaker than our 
standard value <j) w 10~ 2 ; this may occur with a coupling 
(/> oc R proportional to an optical depth, so that the balance 

1 Our referee has kindly pointed out to us that the the new op- 
tical LFs of Richards at al. 2005, preprinted after the submission 
of our MS, agree even better with our predicted excess at faint 
magnitudes. 



condition yields at low a an upturn toward M oc <r 4 oc cr 3 ' 5 . 
The upturn is yet enhanced to M oc <r 3 with a larger scatter, 
when cooling faster than the dynamical time (more likely 
to occur in small dense haloes) offsets the feedback, and 
emulates the condition <f> <C 10~ 2 illustrated by the Fig. lb 
in its left corner, see also CV02. The low a data of Onken 
et al. (2004), and Barth, Greene & Ho (2005) apparently 
indicate these conditions to apply; we expect this may occur 
for a < 70 km/s, an issue that we will develop elsewhere. 

The other imprint left by stochastic but dwindling ac- 
tivity in high density environments is to be found in the 
Eddington ratios that we find to be widely scattered, yet 
declining on average to local values Xe i$ 1/3 for the most 
massive BHs, see Fig. 7. The decline is due primarily to 
the average decrease of the luminosities under the generally 
weakening interactions that tap diminishing gas reservoirs. 
Intrinsic scatter is caused primarily in L and in a milder in- 
tegrated form in M, by the stochastic distribution of orbital 
parameters in the interactions, and by the the dependence 
on a of their effects. Under the control by the feedback the 
overall scatter is actually constrained to under a factor 10 
with some mixing caused by the host dispersion a, as shown 
by the lines in Fig. 7. The result is consistent with the ob- 
servations by McLure & Dunlop (2003) and by Vestergaard 
(2004), considering conceivable sources of additional scatter, 
such as: the initial conditions for z > 2.5 will itself contain 
some dispersion; additional scatter in the data clearly comes 
at low z from observational selection picking up lower value 
of Xe; at high z the current estimates of M grow more un- 
certain based on scaling relations extrapolated from nearby 
AGNs (see Kaspi 2000, Vestergaard 2004). 

Beyond detailed modeling, we stress that the observed 
combination of scattered but declining ratios Xe{z) with 
the tight and steady upper M — a correlation indicates 
supply-limited activation of BHs in stochastic but gener- 
ally dwindling accretion episodes. This is because the later 
and weaker repetitions controlled by the feedback increase 
M and its scatter only moderately, see Fig. 2b; meanwhile 
they decrease the average L and enlarge its variance consid- 
erably, to yield widening scatter to \e(z) superposed to an 
overall declining trend, see Fig. 7. The two figures together 
illustrate how these two trends are made consistent by the 
feedback constraint, and lead us to favour the FCA model 
for 70 < a < 300 km s _1 . 

In sum, a sequence of fueling modes enliven the uniform 
underlying BH paradigm. Out of the several, attendant as- 
trophysical processes our work has focused and linked three 
major components: the dynamical events, that recur with 
decreasing average strength and trigger stochastic but gen- 
erally dwindling accretion episodes; the ensuing depletion 
of the galactic gas reservoirs, that in dense environments 
run out of the supply for accretion and for star formation; 
the constraint imposed to actual accretion by the energy 
feedback from the very source emissions. Our predicted out- 
comes include: flat LFs at higher z > 3, flattening yet at 
fainter luminosities; a steep and sharp M — a correlation for 
the truly supermassive BHs; an upper cutoff to their local 
MD; decline and scatter of the Eddington ratios, widening to 
low z. All features in telling agreement with the developing 
observations. 

To these accretion modes that prevail in dense environ- 
ments and power a first QSO Population, we have added 
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(taking advantage of the linearity of our basic Eqs. 18 and 
25, see §7) in Figs. 6 and 7 also the independent modes 
persisting or standing out for hosts in the field. There we 
have found a second, later AGN population to arise and 
evolve slowly simply because the density-dependent interac- 
tion times r r ~ 1/fijEV given in Eq. (8) are considerably 
longer, around 10 rather than 1 Gyr. Relatedly, high values 
of \e still occur at z ~ 0.5 in some 20% of the faint sources 
and in some 70% of the bright ones; they arise when hosts 
in the field - still gas-rich just because their encounters and 
interactions are late and rare - get involved in strong dy- 
namical events that kindle up relatively small BHs starting 
from under the constrained M — a relation (for related evi- 
dence, see Tanaka 2004). 

We conclude with two straight implications of our pic- 
ture. First, in massive spheroids built up in dense environ- 
ment, the same rapid exhaustion of the gas reservoirs that 
causes strong evolution in the early QSO/AGN population is 
also bound to cause early reddening of the star populations. 
Later and slower field accretion modes, instead, involve blue 
galaxies still rich of gas and actively forming stars. All that 
will be observed (cf. Kauffmann et al. 2003, Hasinger 2004) 
as a bimodal QSO-AGN population, correlated to the bi- 
modal colours of the galaxy population (Dekel & Birnboim 
2004). Second, early small spheroids form small BHs with 
the reduced efficiency discussed in §3 causing the flattened 
and low LF shown in Fig. ; this means a small fraction of BHs 
conspicuous at early z. At low z, instead, many blue galax- 
ies will harbor smaller spheroids and smaller BHs, either 
limited by the feedback constraint or not yet matured; their 
fueling and illumination by field triggering modes - although 
basically driven by the hierarchical formation of the under- 
lying structures - will be perceived as an anti-hierarchical 
development of the activity. 
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Figure 1. la. The bivariate, local distribution N(M, a, z = 0) M 
of the masses M of BHs in hosts with DM velocity dispersion 
cr, when the accretion is unconstrained (the UA model). This is 
obtained after the first, more effective 5 interactions from Eqs. 
(21) and (22), on specifying p(fi \ a) after Eq. (14) (see the end of 
§4). The actual distribution will be smoothed out by the variance 
in the number q of interactions undergone by any given BH. 

lb. The local M — a relation, as obtained on projecting the distri- 
bution in Fig. la onto the M, a plane, see §5 of text. The levels 
shown by grayscale tones are normalized to their maximal value. 
Data points from Ferrarese & Merritt (2000) and from Gebhardt 
et al. (2000). 



© 2005 RAS, MNRAS 000.HHT51 



14 V. Vittorini, F. Shankar, A. Cavaliere 



-1 




2,0 2,2 2,4 2,6 
log ct [km s '] 



Figure 2. Same as Fig. 1, but for feedback-constrained accretion 
(the FCA model). Comparing with Fig. 1, note the tighter relation 
of the averaged M with a. 
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Figure 3. The evolution of the BH mass distribution M N(M,z) 
derived from Eq. (21) for the UA model, and shown at z = 6 
(dashed line), z = 3 (dotted line) and z = (thick solid line). 
The shaded region represents the local estimates given by Yu & 
Tremaine (2002) and by Shankar et al. (2004). 
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Figure 4. Same as Fig. 3, but for the FCA model. Here the thick 
solid line includes satellite captures, while the dot-dashed does 
not. The thin solid line represents the solution of the approximate 
Eq. (19); this yields a close approximation in this case, as expected 
(see text, §5). 
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Figure 5. The optical QSO luminosity function in the form 
L N(L, z) for the UA model; from bottom to top, z = 0.5, 1 and 
2.5. Data from Boyle et al. 2000 and Grazian et al. 2000, marked 
by circles (z = 0.5), triangles (z = 1), and squares (z = 2.5). 
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Figure 6. Same as Fig. 5, for the FCA model. Here we have 
added the LFs at z ~ 0: the thinner solid line represents the con- 
tribution from high density environments (the first QSO/AGN 
population); the dashed lines represent the contributions from low 
density environments (the later AGN population), with the satel- 
lites captures dominating the faint and interactions the brighter 
end. Moreover, the dotted line represents the LF computed with 
the same FCA model at z = 4.5, compared with the observations 
(represented with stars) by Schmidt et al. (1995); Kennefick et 
al. (1995); Fan et al. (2001); Cristiani et al. (2004b). 



0,5-, 




-1,5 - 



-2,0- 

-2,5 I 

0,5 1,0 1,5 2,0 2,5 3,0 

Z 

Figure 7. The evolution of the Eddington ratios Xe computed 
from Eq. (27) in the FCA model, for red host galaxies in high- 
density environments, see §6. The ratios for two specific values 
of a and two probability levels are represented with different line 
styles: the pair of solid lines correspond to <r = 400 km S , and 
the dotted lines to a = 100 km s ; each pair marks the range 
from the zero of Tl q (lower) to the 95% confidence level (upper). 
The redshift z is related to q as given by Eq. (13) and by the 
cosmological t — z relation recalled at the end of §1. The upper 
shading represents the ratios related to blue host galaxies in low- 
density environments, see §7. 
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